% Multisector general equilibrium CalvoPlus model
%
% The state vector is three dimensional. We vectorize it using the
% reshape command. The order of the state variables when they are 
% reshaped is (a, rp, rad)
%
% The structure of this program is very similar to the structure of 
% 'ge1CalvoPlus.m' in the 'One Sector GE' folder. See that file for more 
% detailed documentation.
%
% Jon Steinsson and Emi Nakamura, August 2006
%**************************************************************
function [moments, VC]=geNCalvoPlus_SMM(p0, mu_c, sig_eps_a, rho_a,weight_j );
    
    



% Set Fixed Parameters
%**************************************************************

theta = 4;    % Elasticity of demand
psi = 0;     % One over the Frisch elasticity of labor supply
gamma = 1;   % coefficient of relative risk aversion
FlexSSL = 1/3; % Steady State labor under flexible prices
beta = 0.96^(1/12);   % Discount factor

% Parameters of the process followed by nominal aggregate demand
% log(S_t) - log(S_t-1) = mu + nu_t
% where nu_t ~ N(0,sigma_eta^2), i.i.d.
mu = 0.0028;   %% Inflation trend
sigma_eta = 0.0065; 


mu = 0.00125;   %% Inflation trend
sigma_eta = 0.0028; 




% Set sector specific parameters
%***************************************************************
% NB: Edges of the rad grid need to be extra padded in the multi-sector 
% model because the sectors have different mean rad.
param_agg=[0.00627939500954868 0.0195273437091084];
%param_agg=[0.0501 ;  0.0373 ;   0.0381; 0.6946];
%0.0501    0.0373    0.0381    0.6946
NSectors = 2;
sm = 0.0;
sigma_eps     = [ param_agg(2) sig_eps_a]';
alphaCalvo    = [1-0  1-0 ]';
  % pa           = [p_a    p_a ]';
    rho           = [0.6946    0.6946 ]';
SectorWeights = [ 0.57-weight_j    weight_j]';
%SectorWeights
% rho = 0.4969;
% sigma_eps = 0.0405;
% K = 0.0548/FlexSSY ;
% K2 = K/4000;
% alphaCalvo = 1-0.0349;
% 
% 
% 
% sigma_eps     = [ 0.0397    0.0557      0.0473    0.0405    ]';
% alphaCalvo    = [ 1-0.0903  1-0.0245     1-0.0370  1-0.0349   ]';
% rho           = [0.4942       0.5431     0.5080       0.4969       ]';
% SectorWeights = [ 0.1511    0.1276       0.1146   100   ]';

SectorWeights = SectorWeights/(sum(SectorWeights));
rad_shift =0.01;%0.015
slope = 9; levelshift = 2;
agridsize =40;%40
%%rpgridsize = 3500; %2000;
rpgridsize =400;%300
radgridsize = 30;%25


if abs(sum(SectorWeights)-1)> 10^(-10), error('SectorWeights do not sum to one!'), end

% Maximum number of iterations on G
MaxItOnG = 60;

Ssize = agridsize*rpgridsize*radgridsize;

% Solve for the flexible price steady state 
%****************************************************************

par1 = sm*(theta-1)/(theta-sm*(theta-1));

FlexSSC = FlexSSL*par1^(sm/(1-sm))*(1+par1)^(-1/(1-sm));
FlexSSM = par1*FlexSSC;
FlexSSY = FlexSSC + FlexSSM;

omega = (theta-1)/theta*(1-sm)*(1+par1)*FlexSSL^(-psi-1)*FlexSSC^(1-gamma); % Level parameter on disutility of labor
omega2 = omega*(FlexSSL)^psi; % Marginal disutility of changing prices


K  = [param_agg(2)   mu_c ]'; 

K2 = K/4000;
K = K*FlexSSY;
K2 = K2*FlexSSY;

FlexSSY = log(FlexSSY);
FlexSSL = log(FlexSSL);
FlexSSC = log(FlexSSC);
if sm > 0
    FlexSSM = log(FlexSSM);
end

% % Print Out Parameters
% %****************************************************************
% 
% fprintf('\n')
% fprintf('theta:    %5.5f\n',theta)
% fprintf('psi:      %5.5f\n',psi)
% fprintf('gamma:    %5.5f\n',gamma)
% fprintf('sm:       %5.5f\n',sm)
% fprintf('\n')
% fprintf('Size of State Space: %10.0f\n\n',Ssize)
% 
% fprintf('\n')
% fprintf('K:\n')
% for ii = 1:NSectors
%     fprintf('   %5.4f ',K(ii)/exp(FlexSSY))
% end
% fprintf('\n')
% fprintf('rho:\n')
% for ii = 1:NSectors
%     fprintf('   %5.4f ',rho(ii))
% end
% fprintf('\n')
% fprintf('sigma_eps:\n')
% for ii = 1:NSectors
%     fprintf('   %5.4f ',sigma_eps(ii))
% end
% fprintf('\n\n')

% Set gridsizes and grid max and min
%****************************************************************

MaxUncondVar = max(((sigma_eps.^2)./(1-rho.^2)).^(1/2));
MinUncondVar = min(((sigma_eps.^2)./(1-rho.^2)).^(1/2));
MinUncondVar = max(MinUncondVar,sigma_eta);

trunc_eps = 2.5;
amax = trunc_eps*MaxUncondVar;
amin = -trunc_eps*MaxUncondVar;
ainc = (amax-amin)/(agridsize-1) - mod((amax-amin)/(agridsize-1),10^(-5));
amin = amin - mod(amin,10^(-5));
amax = amin + (agridsize-1)*ainc;
ainc
fprintf('Gridpoints for "a" off of which the least volatile sector is simulated:\n')
fprintf('2*trunc_eps*MinUncondVar/ainc = %5.3f\n\n',...
    2*trunc_eps*MinUncondVar/ainc)
% if 2*trunc_eps*MinUncondVar/ainc < 20
%     error('2*trunc_eps*MinUncondVar/ainc < 20: Results unreliable!',...
%         2*trunc_eps*MinUncondVar/ainc) %#ok<CTPCT>
% end

% rp denotes the real price
% np denotes the nominal price
% The rp grid is a grid for the log of the real price
pimult = 250; % Determines sensitivity of pigrid to mu
rpmax = -amin*(1+pimult*abs(mu));
rpmin = -amax*(1+pimult*abs(mu));
if rpmax < 0.15
    rpmax = 0.15;
    rpmin = -0.15;
end
% %ajout ERwan 18 avril 2017
%  rpmax = 0.3;
%  rpmin = -0.3;

rpinc = (rpmax-rpmin)/(rpgridsize-1) - mod((rpmax-rpmin)/(rpgridsize-1),10^(-5));
rpmin = rpmin - mod(rpmin,10^(-5));
rpmax = rpmin + (rpgridsize-1)*rpinc;

% if rpinc > 0.0075
%     error('rpinc > 0.0075: Results unreliable')
% end

radmin = FlexSSC - radgridsize/2*rpinc + rad_shift;

% Create the real price and idiosyncratic productivity vectors
%**************************************************************

a = amin:ainc:amax; a = a';

rp = zeros(rpgridsize,1);
rp(1) = rpmin;
for ii = 2:rpgridsize
    rp(ii) = rp(ii-1) + rpinc;
end

rad = zeros(radgridsize,1);
rad(1) = radmin;
for ii = 2:radgridsize
    rad(ii) = rad(ii-1) + rpinc;
end

% Create Probability Distribution for the idiosyncratic shock
%**************************************************************
% Proba is a matrix that gives transition probabilities of going
% from state a to state a', i.e. Proba(i,j) = Prob(a' = j | a = i)
% 
% The procedure of Tauchen (1986) is used to approximate the 
% AR(1) process for log(A_it)

for ii = 1:NSectors
   % Proba(ii).Proba = CreateProba_KR(rho(ii),sigma_eps(ii),a, pa(ii));
    Proba(ii).Proba = CreateProba(rho(ii),sigma_eps(ii),a);
end

% Create probability distribution for nominal aggregate demand
%***************************************************************

ProbNgr = CreateProbNgr(mu,sigma_eta,rad);

% Guess a matrix G which gives the log inflation rate given S_t/P_t-1
%**************************************************************

G = zeros(radgridsize,1);
%G = (rp(2)-rp(1))*ones(radgridsize,1);
G(1,1) = -(radgridsize/(2*slope)-mod(radgridsize/(2*slope),1))*(rp(2)-rp(1)) ...
    + levelshift*(rp(2)-rp(1));
for ii = 2:radgridsize
    if mod(ii,slope) == 0
        G(ii,1) = G(ii-1,1) + rp(2) - rp(1);
    else
        G(ii,1) = G(ii-1,1);
    end
end

% Calculate the profit function
%**************************************************************

par1 = psi+1;
par2 = gamma;
par3 = exp(FlexSSC)/exp(FlexSSY);
par4 = exp(FlexSSM)/exp(FlexSSY) - sm;
par5 = 1-sm;

a3d = repmat(a,[1 rpgridsize radgridsize]);

rp3d = repmat(rp,[1 agridsize radgridsize]);
rp3d = permute(rp3d, [2 1 3]);

rad3d = repmat(rad,[1 agridsize rpgridsize]);
rad3d = permute(rad3d,[2 3 1]);

if sm > 0
    l3d = FlexSSL + (par3+par2*par4)/(par5-par1*par4)*(rad3d-FlexSSC);
    m3d = FlexSSM + par1*(l3d - FlexSSL) + par2*(rad3d-FlexSSC);
else
    l3d = FlexSSL + (rad3d-FlexSSC);
    m3d = 0*rad3d;
end

rPi3d = exp(rad3d + m3d).*exp(rp3d).^(1-theta) ...
    - (1-sm)^(sm-1)*sm^(-sm)*omega^(1-sm)*exp(l3d).^(psi*(1-sm)).*exp(rad3d).^(gamma*(1-sm)).*exp(-a3d).*exp(rad3d+m3d).*exp(rp3d).^(-theta);
rPi = reshape(rPi3d,[Ssize 1]);

for ii = 1:NSectors
    menucost(ii).menucost = reshape(omega2*K(ii)*exp(rad3d).^gamma,[Ssize 1]);
end

for ii = 1:NSectors
    menucost2(ii).menucost2 = reshape(omega2*K2(ii)*exp(rad3d).^gamma,[Ssize 1]);
end

clear a3d rp3d rad3d rPi3d l3d m3d

% Calculate the ratio of marginal utility
%**************************************************************
% The (i,j)th element of RMU is the ratio of marginal utility in 
% state j of period t+1 to the marginal utility in state i of
% period t.

RMU = ((1./exp(rad))*(exp(rad)')).^(-gamma);

% Initialize ErV matrix (Expected Value next period)
%**************************************************************

for ii = 1:NSectors
    rV(ii).rV = 1.5*ones(Ssize, 1)*max(0,mean(rPi)/(1-beta));
end

rPi3d = reshape(rPi,[agridsize rpgridsize radgridsize]);
tolV = 0.005*rPi3d(floor(agridsize/2),floor(rpgridsize/2),floor(radgridsize/2))/(1-beta);
clear rPi3d
%tolV = 0.03;

% Initialize stationary distribution Q
%**************************************************************

for ii = 1:NSectors
    Q(ii).Q = ones(Ssize,1)/Ssize;
end

tolQ = Q(1).Q(1)/10;

% Solve for the Equilibrium
%**************************************************************

[G,rV,F,F2,Q] = EqSolverGeNCalvoPlus(rPi,menucost,menucost2,RMU,alphaCalvo,...
    beta,theta,a,rp,rad,Proba,ProbNgr,rV,Q,G,SectorWeights,tolV,tolQ,MaxItOnG,0);

GError = CalcGErrorMultiSectorCalvoPlus(Q,F,F2,G,alphaCalvo,a,rp,rad,theta,...
    Proba,ProbNgr,SectorWeights);

PCStats = RigidityStatsMultiSectorCalvoPlus(Q,F,F2,alphaCalvo,...
    SectorWeights,a,rp,rad);

% 
% fprintf('\n')
% fprintf('GeN CalvoPlus\n\n')
% fprintf('theta:     %5.5f\n',theta)
% fprintf('psi:       %5.5f\n',psi)
% fprintf('gamma:     %5.5f\n',gamma)
% fprintf('sm:        %5.5f\n',sm)
% fprintf('mu:        %5.5f\n',mu)
% fprintf('sigma_eta: %5.5f\n',sigma_eta)
% fprintf('FlexSSC:   %5.5f\n',exp(FlexSSC))
% 
% fprintf('K(i)/FlexSSY:  ')
% for ii = 1:NSectors
%     fprintf('   %5.4f ',K(ii)/exp(FlexSSY))
% end
% fprintf('\n')
% 
% fprintf('K2(i)/FlexSSY: ')
% for ii = 1:NSectors
%     fprintf('   %5.4f ',K2(ii)/exp(FlexSSY))
% end
% fprintf('\n')
% 
% fprintf('1-alphaCalvo:  ')
% for ii = 1:NSectors
%     fprintf('   %5.4f ',1-alphaCalvo(ii))
% end
% fprintf('\n')
% 
% fprintf('rho:           ')
% for ii = 1:NSectors
%     fprintf('   %5.4f ',rho(ii))
% end
% fprintf('\n')
% 
% fprintf('sigma_eps:     ')
% for ii = 1:NSectors
%     fprintf('   %5.4f ',sigma_eps(ii))
% end
% fprintf('\n')
% 
% fprintf('Sector Weights:')
% for ii = 1:NSectors
%     fprintf('   %5.4f ',SectorWeights(ii))
% end
% fprintf('\n')
% 
% fprintf('Overall: \n')
% fprintf('Frequency of Price Change:\n')
% for ii = 1:NSectors
%     fprintf('   %5.4f ',PCStats.freq(ii))
% end
% fprintf('    ||   %5.4f \n',PCStats.avfreq)
% fprintf('Frac Up:\n')
% for ii = 1:NSectors
%     fprintf('   %5.4f ',PCStats.fracup(ii))
% end
% fprintf('    ||   %5.4f \n',PCStats.avfracup)
% fprintf('Size of Price Changes:\n')
% for ii = 1:NSectors
%     fprintf('   %5.4f ',PCStats.wmed(ii))
% end
% fprintf('    ||   %5.4f \n',PCStats.wmed)
% 
% fprintf('Size of Price Changes:\n')
% for ii = 1:NSectors
%     fprintf('   %5.4f ',PCStats.w25(ii))
% end
% fprintf('    ||   %5.4f \n',PCStats.w25)
% fprintf('Size of Price Increases:\n')
% fprintf('Size of Price Changes:\n')
% for ii = 1:NSectors
%     fprintf('   %5.4f ',PCStats.w75(ii))
% end
% fprintf('    ||   %5.4f \n',PCStats.w75)
% fprintf('Size of Price Increases:\n')
% 
% fprintf('    ||   %5.4f \n',PCStats.interq(2))
% fprintf('Interq price changes:\n')
% % for ii = 1:NSectors
% %     fprintf('   %5.4f ',PCStats.sizeup(ii))
% % end
% % fprintf('    ||   %5.4f \n',PCStats.avsizeup)
% % fprintf('Size of Price Decreases:\n')
% % for ii = 1:NSectors
% %     fprintf('   %5.4f ',PCStats.sizedw(ii))
% % end
% % fprintf('    ||   %5.4f \n',PCStats.avsizedw)
% % fprintf('\n')
% % 
% % fprintf('High Cost State\n')
% % fprintf('Frequency of Price Change:\n')
% % for ii = 1:NSectors
% %     fprintf('   %5.4f ',PCStats.freq1(ii))
% % end
% % fprintf('    ||   %5.4f \n',PCStats.avfreq1)
% % fprintf('Frac Up:\n')
% % for ii = 1:NSectors
% %     fprintf('   %5.4f ',PCStats.fracup1(ii))
% % end
% % fprintf('    ||   %5.4f \n',PCStats.avfracup1)
% % fprintf('Size of Price Changes:\n')
% % for ii = 1:NSectors
% %     fprintf('   %5.4f ',PCStats.sizepc1(ii))
% % end
% % fprintf('    ||   %5.4f \n',PCStats.avsizepc1)
% % fprintf('Size of Price Increases:\n')
% % for ii = 1:NSectors
% %     fprintf('   %5.4f ',PCStats.sizeup1(ii))
% % end
% % fprintf('    ||   %5.4f \n',PCStats.avsizeup1)
% % fprintf('Size of Price Decreases:\n')
% % for ii = 1:NSectors
% %     fprintf('   %5.4f ',PCStats.sizedw1(ii))
% % end
% % fprintf('    ||   %5.4f \n',PCStats.avsizedw1)
% % fprintf('\n')
% % 
% % fprintf('Low Cost State\n')
% % fprintf('Frequency of Price Change:\n')
% % for ii = 1:NSectors
% %     fprintf('   %5.4f ',PCStats.freq2(ii))
% % end
% % fprintf('    ||   %5.4f \n',PCStats.avfreq2)
% % fprintf('Frac Up:\n')
% % for ii = 1:NSectors
% %     fprintf('   %5.4f ',PCStats.fracup2(ii))
% % end
% % fprintf('    ||   %5.4f \n',PCStats.avfracup2)
% % fprintf('Size of Price Changes:\n')
% % for ii = 1:NSectors
% %     fprintf('   %5.4f ',PCStats.sizepc2(ii))
% % end
% % fprintf('    ||   %5.4f \n',PCStats.avsizepc2)
% % fprintf('Size of Price Increases:\n')
% % for ii = 1:NSectors
% %     fprintf('   %5.4f ',PCStats.sizeup2(ii))
% % end
% % fprintf('    ||   %5.4f \n',PCStats.avsizeup2)
% % fprintf('Size of Price Decreases:\n')
% % for ii = 1:NSectors
% %     fprintf('   %5.4f ',PCStats.sizedw2(ii))
% % end
% % fprintf('    ||   %5.4f \n',PCStats.avsizedw2)
% % fprintf('\n')
% % 
% % fprintf('Fraction of Price Changes in Low Cost State:\n')
% % for ii = 1:NSectors
% %     fprintf('   %5.4f ',(1-alphaCalvo(ii))*PCStats.freq2(ii)/PCStats.freq(ii))
% % end

% Save Model Input and Output
%**************************************************************
% 
% Inputs
GENModel_med2_GLB.a = a;
GENModel_med2_GLB.rp = rp;
GENModel_med2_GLB.rad = rad;
GENModel_med2_GLB.K = K;
GENModel_med2_GLB.K2 = K2;
GENModel_med2_GLB.beta = beta;
GENModel_med2_GLB.theta = theta;
GENModel_med2_GLB.sm = sm;
GENModel_med2_GLB.gamma = gamma;
GENModel_med2_GLB.omega = omega;
GENModel_med2_GLB.omega2 = omega2;
GENModel_med2_GLB.rho = rho;
GENModel_med2_GLB.sigma_eps = sigma_eps;
GENModel_med2_GLB.mu = mu;
GENModel_med2_GLB.sigma_eta = sigma_eta;
GENModel_med2_GLB.SectorWeights = SectorWeights;
GENModel_med2_GLB.alphaCalvo = alphaCalvo;

% Outputs
GENModel_med2_GLB.G = G;
%GENModel.rV = rV;
GENModel_med2_GLB.F = F;
GENModel_med2_GLB.F2 = F2;
GENModel_med2_GLB.Q = Q;
GENModel_med2_GLB.GError = GError;
GENModel_med2_GLB.freq = PCStats.freq1;
GENModel_med2_GLB.fracup = PCStats.fracup1;
GENModel_med2_GLB.sizepc = PCStats.sizepc1;
GENModel_med2_GLB.sizeup = PCStats.sizeup1;
GENModel_med2_GLB.sizedw = PCStats.sizedw1;
GENModel_med2_GLB.avfreq = PCStats.avfreq1;
GENModel_med2_GLB.avfracup = PCStats.avfracup1;
GENModel_med2_GLB.avsizepc = PCStats.avsizepc1;
GENModel_med2_GLB.avsizeup = PCStats.avsizeup1;
GENModel_med2_GLB.avsizedw = PCStats.avsizedw1;
GENModel_med2_GLB.freq = PCStats.freq2;
GENModel_med2_GLB.fracup = PCStats.fracup2;
GENModel_med2_GLB.sizepc = PCStats.sizepc2;
GENModel_med2_GLB.sizeup = PCStats.sizeup2;
GENModel_med2_GLB.sizedw = PCStats.sizedw2;
GENModel_med2_GLB.avfreq = PCStats.avfreq2;
GENModel_med2_GLB.avfracup = PCStats.avfracup2;
GENModel_med2_GLB.avsizepc = PCStats.avsizepc2;
GENModel_med2_GLB.avsizeup = PCStats.avsizeup2;
GENModel_med2_GLB.avsizedw = PCStats.avsizedw2;
GENModel_med2_GLB.freq = PCStats.freq;
GENModel_med2_GLB.fracup = PCStats.fracup;
GENModel_med2_GLB.sizepc = PCStats.sizepc;
GENModel_med2_GLB.sizeup = PCStats.sizeup;
GENModel_med2_GLB.sizedw = PCStats.sizedw;
GENModel_med2_GLB.avfreq = PCStats.avfreq;
GENModel_med2_GLB.avfracup = PCStats.avfracup;
GENModel_med2_GLB.avsizepc = PCStats.avsizepc;
GENModel_med2_GLB.avsizeup = PCStats.avsizeup;
GENModel_med2_GLB.avsizedw = PCStats.avsizedw;
GENModel_med2_GLB.avwmed=PCStats.avwmed;
GENModel_med2_GLB.interq=PCStats.avinterq;
GENModel_med2_GLB.avkur=PCStats.avkur;
save('GENModel_med2_GLB','GENModel_med2_GLB')

%toc
 

Qrad = zeros(radgridsize,1);
for ii = 1:NSectors
    Qtemp = reshape(Q(ii).Q,[agridsize*rpgridsize radgridsize]);
    Qrad = Qrad + SectorWeights(ii)*sum(Qtemp)';
end

Erad = rad'*Qrad;
VARrad = Qrad'*((rad-Erad).^2);

% Print Results
%**************************************************************************

fprintf('St.Dev of output = %5.5f (x10^-2)\n',(VARrad^(1/2))*100)
fprintf('Variance of output = %5.9f (x10^-4)\n\n',VARrad*10000)
 


 f2=PCStats.freq(2);
   avfracup2=PCStats.fracup(2);
  med2=PCStats.wmed(2);
  q12=PCStats.w25(2);
  q32=PCStats.w75(2);
  interq2=PCStats.interq(2);
  kur2=PCStats.kur(2);
    share_Calvo=(1-alphaCalvo(2))*PCStats.freq2(2)/PCStats.freq(2);
  mu_p=mu_c*3/4;
  MC_paid=f2*(1-share_Calvo)*mu_p*100;
  moments=[f2  avfracup2 med2 interq2 kur2 q12 med2 q32 share_Calvo mu_p MC_paid];
    VC=VARrad*100;

end
